The following sections will cover how to pull out common diagnostics from VAST model outputs and some data visualization development for cross boundary metrics for US and Canadian population densities.
EDIT: These outputs will be working from SDMTMB results and not VAST. This means I will also be demoing new post-processing code.
# Load packageslibrary(units)
udunits database from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/units/share/udunits/udunits2.xml
library(here)
here() starts at /Users/adamkemberling/Documents/Repositories/Coca_Shiny
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(gmRi)library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
Registered S3 method overwritten by 'ggside':
method from
+.gg ggplot2
# Paths to Box Assets mills_path <-cs_path(box_group ="mills")coca_path <-str_c(mills_path, "Projects/COCA19_Projections/")projections_path <-paste0(coca_path, "projections/")mods_path <-paste0(coca_path, "mod_fits/")# Hexagon gridhex_grid <-read_sf(here::here("COCA_SDM_app_dev/dev/scratch_data", "hex_grid.geojson"))# Load the shapefilesdfo_bounds <-read_sf(here::here("local_data/Regions_for_CRSBND/DFO.shp"))nmfs_bounds <-read_sf(here::here("local_data/Regions_for_CRSBND/NMFS.shp"))land_sf <-read_sf(here::here("COCA_SDM_app_dev/dev/scratch_data", "nw_atlantic_countries_crs32619.geojson"))# Load the Hague Lineshague_sf <-read_sf(here::here("COCA_SDM_app_dev/dev/scratch_data", "hagueline_crs32619.geojson"))# Need to get an area per region...number of cells per regionall_regs <-bind_rows(dfo_bounds, nmfs_bounds) |>mutate(area =st_area(geometry)) %>%rename(jurisdiction = Region)all_regs$area<-set_units(all_regs$area, km^2)
# Plotting map themetheme_map <-function(fontfam ="Avenir", guides = T, ...){list(# Theme options, with ellipse to add moretheme(# Font across all texttext =element_text(family ="Avenir"),# Titles + Textplot.title =element_text(hjust =0, face ="bold", size =20),plot.subtitle =element_text(size =18),legend.title =element_text(size =16, lineheight =1.75),legend.text =element_text(size =14), legend.spacing.y =unit(1.75, "lines"),# Grids and Axespanel.background =element_blank(), panel.border =element_rect(color ="black", fill ="transparent"), panel.grid.major =element_line(color ="gray80"),axis.title.x =element_blank(),axis.title.y =element_blank(),axis.ticks=element_blank(),plot.margin =margin(t =1, r =1, b =1, l =1, unit ="pt"),#legend.position = c(.725, .125), legend.background =element_rect(color ="transparent", fill ="white", linewidth =0.25),# Facetsstrip.background =element_rect(fill ="#00736D", color ="white"),strip.text =element_text(color ="white", face ="bold",size =12,family = fontfam),legend.position ="bottom",# Use ellipses for tweaks on the fly: ...))}
# This was found in Create_Mesh_from_Knots.Rsf_meshify <-function(input_df, coords =c("Lon", "Lat"), length_km =25, in_crs =4326, trans_crs =32619, square = T){# Make the dataframe an sf class using coordinates in_sf <-st_as_sf(input_df, coords = coords, crs = in_crs, remove = F) %>%# Transform it to a crs that is projected in metersst_transform(crs = trans_crs)# If we are getting gaps we can buffer here, or up the mesh size# Use that data to define a grid with dimensions of length_km*length_km sf_grid <-st_make_grid(x = in_sf,cellsize =c(length_km*1000, length_km*1000), what ="polygons", square = square) %>%# Make the grid an sf classst_as_sf() # Use the original data to trim it so its just cells that overlap the points sf_out <- sf_grid %>%st_filter(in_sf, .predicate = st_contains) %>%st_as_sf() # Join the clipped grid to the dataset sf_out <-st_join(sf_out, in_sf, join = st_intersects)# Return the results return(sf_out)}
# Loading the VAST Modelsmod_list <-setNames(list.files(mods_path, pattern =".rds", full.names = T),str_remove(list.files(mods_path, pattern =".rds"), ".rds"))# Load a random modelmod_test <-readRDS(mod_list[[2]])# OR, use haddock / spiny dogfish as an example had_vast <-read_rds(str_c(mods_path, "Haddock_full_fitted_vast.rds"))sd_vast <-read_rds(str_c(mods_path, "SpinyDogfish_full_fitted_vast.rds"))# Get some details on the range parametersdens_dep_range <- sd_vast$Report$Range_raw1 # density dependent rangedens_ind_range <- sd_vast$Report$Range_raw2 # density independent range
VAST Model Diagnostics
Each VAST model is estimating a common suite of terms which we can check and investigate across species. Every model will have the following three statistical relationships: * Long-term spatial autocorrelation random effect. Adjusts presence/absence and density estimates by location * Seasonal adjustments
# This was an unnecessary exercise, knots are not where distribution data is located#------- Mapping with the hex_grid# Make a hex mesh from the know coordinatesknot_hex_grid <-sf_meshify(input_df =distinct(spat_grid, Lon, Lat), square = F, length_km =35)
p1_map <- spat_grid %>%st_drop_geometry() %>% dplyr::filter(omega =="Presence_Absence") %>%# Join the hexagon meshleft_join(knot_hex_grid, by =join_by(Lon, Lat)) %>%st_as_sf() %>%ggplot() +geom_sf(aes(fill = omega_val)) +facet_wrap(~omega) +scale_fill_distiller(palette ="RdBu") +theme_dark() +theme(legend.title.position ="top", legend.title =element_text(hjust =0.5),) +labs(fill ="Presence/Absence Random Effect")p2_map <- spat_grid %>%st_drop_geometry() %>%filter(omega =="Density") %>%# Join the hexagon meshleft_join(knot_hex_grid, by =join_by(Lon, Lat)) %>%st_as_sf() %>%ggplot() +geom_sf(aes(fill = omega_val)) +facet_wrap(~omega) +scale_fill_distiller(palette ="RdBu") +theme_dark() +theme(legend.title.position ="top", legend.title =element_text(hjust =0.5),) +labs(fill ="Abundance Density Random Effect")p1_map / p2_map +plot_annotation(title ="Spiny Dogfish Spatial Random Effects")
These models contain three environmental covariates: surface temperature, bottom temperature, and depth. Each of these covariates will have an associated preference curve showing how species cpue changes along a continuous range of these values.
test_prefs <-filter(pref_all, comname =="spiny dogfish")test_prefs %>%ggplot() +geom_area(aes(val_actual, y = fit_exp, group = comname), alpha =0.4, fill =gmri_cols("blue economy teal")) +geom_line(aes(val_actual, fit_exp, group = comname), linewidth =1) +facet_wrap(.~variable, scales ="free", ncol =2) +scale_color_gmri() +guides(linetype =guide_legend(title.position ="top", title.hjust =0.5),color =guide_legend(title.position ="top", title.hjust =0.5)) +theme_minimal() +labs(x ="Covariate Value", title ="Spiny Dogfish Preferences", color ="Average Regional Condition",linetype ="SSP Scenario",y ="Predicted kg/km2")
The code below digs into how Andrew Actually accomplished extracting these details from the VAST models. Its a bit much for this workup…
VAST Post-Processing of Density Estimates
# Load the Projections / Model Outputs - Raw# I hate nested dataframes, so I'm just going to use a list at each step:read_rds_func <-function(file_name) { out <-readRDS(paste0(file_name))return(out)}# Loads all Files in the folder that end with _mean.rdsfpaths <-list.files(projections_path, pattern ="_mean.rds", full.names =TRUE)fnames <-str_remove(list.files(projections_path, pattern ="_mean.rds", full.names =FALSE, ), ".rds")proj_files <-setNames(fpaths, fnames)# Read each dataset into a list or a testerprojection_list <-map(proj_files, read_rds_func)projection_test <-read_rds_func(proj_files["SpinyDogfish_full_CMIP6_SSP5_85_mean"])#### 2. Isolate the Biomass Density Information ##### For now, we are mostly going to be using the results in the "Dens" object to make our maps # and then also for cropping with community footprints to get our summaries of change. # Let's pull out the density results as a new column# Isolate just the density informationdensity_estimates <- projection_list %>%map(~pluck(.x, "Dens"))density_test <- density_estimates[["SpinyDogfish_full_CMIP6_SSP5_85_mean"]] %>%mutate(season =case_match(month(Time), 3~"Spring",7~"Summer", 10~"Fall"),season =factor(season, levels =c("Spring", "Summer", "Fall"))) %>%ungroup()
Mapping VAST Projections
# Make a hex mesh from the know coordinatespredictions_hex_grid <-sf_meshify(input_df =distinct(density_test, Lon, Lat), square = F, length_km =35) # Plot thoselibrary(scales)
Attaching package: 'scales'
The following object is masked from 'package:purrr':
discard
The following object is masked from 'package:readr':
col_factor
density_test %>%filter(year(Time) ==2025) %>%left_join(predictions_hex_grid) %>%st_as_sf() %>%#st_as_sf(coords = c("Lon", "Lat"), crs = 4326, remove = F) %>% ggplot() +geom_sf(aes(fill = Prob_0.5)) +scale_fill_distiller(palette ="RdBu", transform ="log10", labels =label_log()) +facet_wrap(~season, ncol =2) +theme_dark() +labs(title ="Spiny Dogfish 2025 Projections: SSP5-8.5 Mean Estimate",subtitle ="Vast Density Estimates use the projection grid, not knots")
Joining with `by = join_by(Lat, Lon)`
SDMTMB Results
# Andrew sent three models related to lobsterlob_juve <-read_rds(here::here("local_data","Example Projection Files", "juve_projected_biomass.rds"))lob_adult <-read_rds(here::here("local_data","Example Projection Files", "adult_projected_biomass.rds"))lob_pred <-read_rds(here::here("local_data","Example Projection Files", "pred_projected_biomass.rds"))lob_projections <-bind_rows(list("Juvenile Lobsters"= lob_juve,"Adult Lobsters"= lob_adult#,#"Lobster Predators" = lob_pred), .id ="model_group") %>%mutate(season =factor(season, levels =c("Spring", "Summer", "Fall")),model_group =factor(model_group, levels =c("Juvenile Lobsters", "Adult Lobsters", "Lobster Predators")) )
Proportion in Each Region
For cross-boundary we are interested in differences between US and Canadian study areas.
# What is the proportion of the predicted biomass in each region# # Do the overlays to assign which region it belongs too# dfo_bounds# nmfs_bounds# Get unique Coordinatessdm_locations <-distinct(lob_projections, longitude, latitude) %>%st_as_sf(coords =c("longitude", "latitude"), crs =4326, remove = F)# Make a mesh from the sdm locationssdm_hex_grid <-sf_meshify(input_df =distinct(lob_projections, Lon = longitude, Lat = latitude), square = F, length_km =35)# Label which area the points fall withindfo_locations <-st_join(sdm_locations, dfo_bounds, join = st_within) %>%filter(is.na(Region) ==FALSE)nmfs_locations <-st_join(sdm_locations, nmfs_bounds, join = st_within) %>%filter(is.na(Region) ==FALSE)region_labs <-bind_rows(dfo_locations, nmfs_locations) %>%st_drop_geometry() %>%rename(jurisdiction = Region)# Add these back inlob_projections <-left_join(lob_projections, region_labs)
Joining with `by = join_by(longitude, latitude)`
Baseline Period Maps
# Get the average cpue over the most recent 20 yearslob_baselines <- lob_projections %>%filter(year %in%c(2004:2023)) %>%group_by(jurisdiction, model_group, season, longitude, latitude) %>%summarise(baseline_biomass_mean =exp(mean(proj_biomass_mean, na.rm = T)),.groups ="drop")# Map the Baseline Period CPUElob_baselines %>%rename(Lon = longitude, Lat = latitude) %>%left_join(sdm_hex_grid) %>%st_as_sf() %>%ggplot() +geom_sf(aes(fill = baseline_biomass_mean), alpha =0.75) +geom_sf(data = land_sf, color ="gray95", fill ="gray40", linewidth =0.15) +geom_sf(data = hague_sf, color ="black", linewidth =1, linetype =1) + rcartocolor::scale_fill_carto_c(palette ="RedOr", transform ="log10",labels =label_number(accuracy =0.01)) +facet_grid(model_group~season) +theme_map() +theme(legend.position ="bottom", legend.title.position ="top") +guides(fill =guide_colorbar(barwidth =unit(10, "cm"))) +coord_sf(xlim =c(-182500, 1550000), ylim =c(3875000, 5370000) , expand = F, crs =32619) +labs(fill ="Baseline Biomass CPUE kg/km2")
Joining with `by = join_by(Lon, Lat)`
Jurisdictional Timeseries
####### Biomass time series###### Average density within each region per time stepres_summ <- lob_projections |>group_by(season, year, Date, jurisdiction) |>summarize("mean_biomass"=mean(proj_biomass_mean)) %>%filter(is.na(jurisdiction) ==FALSE)
`summarise()` has grouped output by 'season', 'year', 'Date'. You can override
using the `.groups` argument.
The catch/effort change from one period to another could be compared as either percent change or log ratio (log fold change). Here is what those could look like
For these we will have the outputs from two models, so I am going to reshape the data and treat the juvenile lobster densities as one model and the adult lobster densities as another.
x axis = biomass density change y = year with jurisdiction offset
# Average density within each region per time stepperiod_summ <- lob_projections |>mutate(period =case_when( year %in%c(2004:2023) ~"base",#year %in% c(2010:2023) ~ "base", year %in%c(2046:2055) ~"mid", year %in%c(2091:2100) ~"end",TRUE~"drop" )) %>%filter(period !="drop") %>%group_by(model_group, season, period, jurisdiction) |>summarize("mean_biomass"=mean(proj_biomass_mean),.groups ="drop") %>%filter(is.na(jurisdiction) ==FALSE) |># Join, calculate total arealeft_join(all_regs, by =c("jurisdiction")) |>mutate("Biomass_Total"=drop_units(exp(mean_biomass)*area)) %>%select(-c(geometry, area))# Reshapeperiod_differences <- period_summ %>%select(-mean_biomass) %>%pivot_wider(names_from = period, values_from = Biomass_Total) %>%mutate(mid_change = mid - base,end_change = end - base,mid_change_percent = mid_change/base,end_change_percent = end_change/base,#season = fct_rev(season),scenario ="SSP5: 8.5")
# Build a Plotggplot(period_differences) +geom_vline(xintercept =0) +geom_point(aes(mid_change_percent, y = season, color = model_group),size =2) +scale_x_continuous(labels =label_percent()) +scale_color_gmri() +theme_classic() +facet_grid(scenario~.) +labs(x ="Change in Biomass",y ="Season")
#---------# Initial Thoughts, I think I may want actuall Biomass?# but even if I don't I think I want the different periods on their own lines# lets try chicklets# Do some major reshapingchicklet_prep <- period_differences %>%select(-c(base, mid, mid_change, end, end_change)) %>%pivot_longer(cols =ends_with("percent"), names_to ="period", values_to ="percent_change") %>%mutate(model_group =if_else( model_group =="Juvenile Lobsters", "env_only", "vast"),bar_y =if_else(str_detect(period, "mid"), 1.1, 1),period_label =if_else(str_detect(period, "mid"), "Mid-Century", "Turn of the Century"),period_alpha =if_else(str_detect(period, "mid"), 1, 0.7)) %>%pivot_wider(names_from = model_group, values_from ="percent_change") %>%mutate(label_x = (env_only + vast)/2)#chicky_points <- period_differences %>% # Plot thatbar_height <-0.025ggplot(chicklet_prep) +geom_vline(xintercept =0, linewidth =0.8, color ="gray60") + ggchicklet:::geom_rrect(aes(xmin = vast,xmax = env_only, ymin = bar_y - bar_height, ymax = bar_y + bar_height, fill = period_label,alpha =I(period_alpha)),r =unit(0.45, 'npc')) +facet_nested( jurisdiction + season ~ scenario, nest_line =TRUE) +scale_x_continuous(limits =c(-1,1), labels =label_percent()) +scale_y_continuous(expand =expansion(mult =c(0.75,0.75))) +scale_fill_gmri() +theme_minimal() +theme(text =element_text(family ="Avenir"),legend.position ="top",legend.title.position ="left",axis.text.y =element_blank(),legend.title =element_text(size =12, face ="bold"),panel.grid.major.y =element_blank(),strip.text =element_text(family ="Avenir", face ="bold", size =12)) +labs(x ="Biomass Change",fill ="Projection Period", y ="",title ="Projected Regional Biomass Changes",subtitle ="Percent difference from 2004-2023 Baseline Period")
# Do some major reshapingchicklet_prep <- period_differences %>%select(-c(base, mid, mid_change, end, end_change)) %>%pivot_longer(cols =ends_with("percent"), names_to ="period", values_to ="percent_change") %>%mutate(model_group =if_else( model_group =="Juvenile Lobsters", "env_only", "vast"),bar_y =if_else(jurisdiction =="DFO", 1, 1.1),period_label =if_else(str_detect(period, "mid"), "Mid-Century", "Turn of the Century")) %>%pivot_wider(names_from = model_group, values_from ="percent_change") %>%mutate(label_x = (env_only + vast)/2)# Plot thatbar_height <-0.025ggplot(chicklet_prep) +geom_vline(xintercept =0, linewidth =0.8, color ="gray60") + ggchicklet:::geom_rrect(aes(xmin = vast,xmax = env_only, ymin = bar_y - bar_height, ymax = bar_y + bar_height, fill = jurisdiction),r =unit(0.45, 'npc')) +geom_point(data = chicklet_prep,aes(x = env_only, y = bar_y, shape ="Environment\nOnly Model"), size =3, fill ="white") +geom_point(data = chicklet_prep,aes(x = vast, y = bar_y, shape ="Spatio-Temporal\nModel"), size =3, fill ="white") +facet_nested( period_label + season ~ scenario, nest_line =TRUE) +scale_shape_manual(values =c(21, 24)) +scale_x_continuous(limits =c(-1,1), labels =label_percent()) +scale_y_continuous(expand =expansion(mult =c(0.75,0.75))) +scale_fill_gmri() +theme_minimal() +theme(text =element_text(family ="Avenir"),legend.position ="bottom",legend.title.position ="left",axis.text.y =element_blank(),legend.title =element_text(size =12, face ="bold"),panel.grid.major.y =element_blank(),strip.text =element_text(family ="Avenir", face ="bold", size =12),legend.box ="Vertical") +labs(x ="Biomass Change",fill ="National Jurisdiction:", y ="",title ="Projected Regional Biomass Changes",subtitle ="Percent difference from 2004-2023 Baseline Period",shape ="Distribution Model Type:")
# Can we make a legend that explains what the bars mean?
Center of Gravity Figure
For this section I want to use density curves in the margins to capture decadal scale movements in the centers of gravity
# Center of gravity as weighted meanCOGrav<-function(df) { cog_lat<-weighted.mean( df[,"latitude"], w =exp(df[, "proj_biomass_mean"])) cog_lon<-weighted.mean( df[,"longitude"], w =exp(df[, "proj_biomass_mean"])) cog_out<-data.frame("COGx"= cog_lon, "COGy"= cog_lat)return(cog_out)}# Estimate center of gravityres_cog <- lob_projections |>group_by(model_group, season, year) |>nest() %>%mutate(ID =row_number(),COG =map(data, ~COGrav(.x)) ) |># dplyr::select(year, season, COG) |>unnest_wider(COG) %>%mutate(season =factor(season, levels =c("Spring", "Summer", "Fall")),Decade = year - year %%10,Decade =factor(Decade, levels =sort(unique(Decade))))
# Does not work with sf, Just use year-roundlibrary(ggside)# Estimate center of gravityres_cog_annual <- lob_projections |>group_by(model_group, year) |>nest() %>%mutate(ID =row_number(),COG =map(data, ~COGrav(.x)) ) |># dplyr::select(year, season, COG) |>unnest_wider(COG) %>%mutate(Decade = year - year %%10,Decade =factor(Decade, levels =sort(unique(Decade))))# ggplot(res_cog_annual, aes(COGx, COGy, color = Decade)) +geom_point(aes(color = Decade), alpha =0.9) +geom_xsidedensity(aes(y =after_stat(density), fill = Decade), position ="dodge", alpha =0.4) +geom_ysidedensity(aes(x =after_stat(density), fill = Decade), position ="dodge", alpha =0.4) +facet_grid( model_group~.,labeller =labeller(model_group =label_wrap_gen(width =8))) +theme_classic() +theme(text =element_text(family ="Avenir"),legend.position ="bottom",legend.title.position ="left",axis.text.y =element_blank(),legend.title =element_text(size =12, face ="bold"),panel.grid.major.y =element_blank(),strip.text =element_text(family ="Avenir", face ="bold", size =12)) +theme(panel.grid.major =element_blank()) +guides(color =guide_legend(override.aes =list(shape =15, size =5),nrow =2)) +theme(ggside.panel.scale = .3) +labs(title ="Center of Biomass Shifts",subtitle ="Lat + Lon Components Displayed as Decadal Distributions",y ="Latitude",x ="Longitude")
Warning: Groups with fewer than two data points have been dropped.
Groups with fewer than two data points have been dropped.
Groups with fewer than two data points have been dropped.
Groups with fewer than two data points have been dropped.
Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`
Width not defined
ℹ Set with `position_dodge(width = ...)`
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
# Could I put a little flag pole or tick mark where the decadal centers are? Yesdecadal_cogs <- res_cog_annual %>%group_by(Decade, model_group) %>%summarise(across(starts_with("COG"), mean))
`summarise()` has grouped output by 'Decade'. You can override using the
`.groups` argument.
new_england <- land_sf %>%st_transform(st_crs(4326))ggplot(res_cog_annual) +geom_point(data = res_cog_annual,aes(COGx, COGy, color = Decade), alpha =0.9) +geom_xsidevline(data = decadal_cogs,aes(xintercept = COGx, color = Decade), linewidth =1, alpha =0.9) +geom_ysidehline(data = decadal_cogs,aes(yintercept = COGy, color = Decade), linewidth =1, alpha =0.9) +# # Turns off ticks for side axes, not working with sf# scale_xsidey_continuous(# breaks = NULL, labels = "", expand = expansion(c(0,.1))) +# scale_ysidex_continuous(# breaks = NULL, labels = "", expand = expansion(c(0,.1))) +# What about land now?geom_sf(data = new_england, color ="dark gray", lwd =0.2, na.rm =TRUE) +facet_wrap(~model_group, ncol =2,labeller =labeller(model_group =label_wrap_gen(width =8))) +# Use scale_x_continuous to zoom, coord_sf is an issuescale_color_carto_d(palette ="Vivid") +scale_x_continuous(limits =range(res_cog_annual$COGx) +c(-2,2)) +scale_y_continuous(limits =range(res_cog_annual$COGy) +c(-1,1)) +guides(color =guide_legend(override.aes =list(shape =15, size =5),nrow =2)) +theme_classic() +theme(text =element_text(family ="Avenir"),legend.position ="bottom",legend.title.position ="left",axis.text.x =element_blank(),axis.text.y =element_blank(),axis.ticks.x =element_blank(),axis.ticks.y =element_blank(),legend.title =element_text(size =12, face ="bold"),panel.grid.major.y =element_blank(),panel.grid.major =element_blank(),strip.text =element_text(family ="Avenir", face ="bold", size =12)) +theme(ggside.panel.scale = .05) +labs(title ="Center of Biomass Shifts",subtitle ="Decadal Centerpoints Displayed in Margins",y ="Latitude",x ="Longitude")
decadal_cog_ranges <- res_cog_annual %>%group_by(Decade, model_group) %>%summarise(across(starts_with("COG"), .fns =list("min"= min, "max"= max)),.groups ="drop") %>%mutate(decade_y =as.numeric(Decade))ggplot(res_cog_annual) +geom_point(data = res_cog_annual,aes(COGx, COGy, color = Decade), alpha =0.9) +geom_xsidesegment(data = decadal_cog_ranges,aes(x = COGx_min, xend = COGx_max, y = decade_y, yend = decade_y, color = Decade), linewidth =1, alpha =0.9) +geom_ysidesegment(data = decadal_cog_ranges,aes(x = decade_y, xend = decade_y, y = COGy_min, yend = COGy_max, color = Decade), linewidth =1, alpha =0.9) +# # Turns off ticks for side axes, not working with sf# scale_xsidey_continuous(# breaks = NULL, labels = "", expand = expansion(c(0,.1))) +# scale_ysidex_continuous(# breaks = NULL, labels = "", expand = expansion(c(0,.1))) +# What about land now?geom_sf(data = new_england, color ="dark gray", lwd =0.2, na.rm =TRUE) +facet_wrap(~model_group, ncol =2,labeller =labeller(model_group =label_wrap_gen(width =8))) +# Use scale_x_continuous to zoom, coord_sf is an issue#scale_color_carto_d(palette = "Vivid") +scale_color_carto_d(palette ="Mint", direction =-1) +scale_x_continuous(limits =range(res_cog_annual$COGx) +c(-2,2)) +scale_y_continuous(limits =range(res_cog_annual$COGy) +c(-1,1)) +guides(color =guide_legend(override.aes =list(shape =15, size =5),nrow =2)) +theme_classic() +theme(text =element_text(family ="Avenir"),legend.position ="bottom",legend.title.position ="left",axis.text.x =element_blank(),axis.text.y =element_blank(),axis.ticks.x =element_blank(),axis.ticks.y =element_blank(),legend.title =element_text(size =12, face ="bold"),panel.grid.major.y =element_blank(),panel.grid.major =element_blank(),strip.text =element_text(family ="Avenir", face ="bold", size =12)) +theme(ggside.panel.scale = .1) +labs(title ="Center of Biomass Shifts",subtitle ="Decadal Centerpoints Displayed in Margins",y ="Latitude",x ="Longitude")
